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Radio observations of neutron star binary pulsar systems have constrained strongly the masses of eight neutron stars. Assuming 
neutron star masses are uniformly distributed between lower and upper bounds mi and m u , the observations determine with 
95% confidence that 1.01 < mi/Mg < 1.34 and 1.43 < m„/M < 1.64. These limits give observational support to neutron 
star formation scenarios that suggest that masses should fall predominantly in the range 1.3 < m/M < 1.6, and will also be 
important in the interpretation of binary inspiral observations by the Laser Interferometer Gravitational-wave Observatory. 
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a. Introduction and Motivation Radio observations 
of four neutron star - neutron star (ns-ns) binary pul- 
sar systems have constrained the masses of eight neutron 
stars, and all the masses lie close to 1.35 Mq. This coin- 
cidence suggests that natural formation mechanisms re- 
strict the range of ns masses more than limitations due 
to the nuclear and super-nuclear equation of state. This 
suggestion also arises in theoretical studies of core col- 
lapse supernovae [1]. Here I explore the statistical signif- 
icance of this coincidence, with the goal of determining 
nature's limits on ns masses: in particular, modeling the 
ns mass distribution as uniform between upper and lower 
bounds m u and m;, I use the observations to determine 
the joint probability distribution of m u and m;. 

While noted before [2], the coincidence in measured ns 
masses has never been subjected to a statistical analy- 
sis that treats the observations jointly and accounts for 
the statistical uncertainties in the mass determinations. 
Here I use Bayesian statistical techniques to combine the 
separate observations and determine the probability dis- 
tribution of the mass bounds. 

A second motivation for this study is the anticipated 
observation of binary ns inspiral by the Laser Interferom- 
eter Gravitational- wave Observatory (LIGO) [3]. LIGO 
will be sensitive to ns masses in inspiraling binaries [4,5], 
and interpretation of its observations will require an ac- 
curate assessment of our knowledge of the mass distribu- 
tion in ns-ns binaries [6]. This work is meant to begin 
that assessment. 

b. Summary of relevant observations Observations of 
binary pulsars 1913+16 and 1534+12 have determined 
for each system the total mass M and the pulsar com- 
panion mass m c [7,8]. In both cases the companion is 
believed to be a ns. Similarly, observations of binary 
pulsars 2127+11C and 2303+46 have determined for each 
M and the mass function / [9,2]. Assuming a probability 
density P(cosi) = 1/2 for the orbital inclination angle i 
limits the mass of the pulsar and its companion, which 
(for these two systems) is believed to be a ns. Here I 
assume that the measured /, M , and m c are distributed 



TABLE I. The values adopted for the total mass M, the 
companion mass m c , and the standard error of each (<tm and 
<Tm c ) of PSR1913+16 and PSR1534+12 [7,8]. 



System M/M 


<tm/M m c /M 


<Tm c /M 


1913+16 2.82827 
1534+12 2.679 


4 x 10" 5 1.442 
0.003 1.36 


0.003 
0.03 


TABLE II. The values adopted for the mass function /, 
the total mass M, and the standard error of each (<r/ and 
<jm) of PSRs 2127+11C and 2303+46 [9,2,17]. 


System //M 


aj/M @ M/M 


<7 M /M 


PSR2127+11C 0.15285 
PSR2303+46 0.246287 


1.8 X 10 -4 2.706 
6.7 X 10" 6 2.57 


3.6 X 10 -3 
0.08 



independently and normally about the actual /, M , and 
fh c , where the distribution variances are given by the re- 
ported 1 a uncertainties in /, M , and m c : 



P(x\x, T) = exp 



/V27T(T X 



(1) 



where x is one of /, M , or m c , and X represents other, 
unenumerated assumptions that characterize the obser- 
vations. Tables I and II give the values adopted here for 
/, m c , M, a j, <7 mc and ctm of these four systems. 

The relationship between /, M , m c , /, M , and m c is 
not as simple as equation (1); however, for small vari- 
ances the Gaussian approximation is a good one, and 
in the absence of a detailed description of the fitting 
procedure that determines /, m c , and M from the ob- 
servations the assumption of statistical independence is 
reasonable. Additionally, numerical investigations show 
that the final 95% confidence intervals for m; and m u are 
largely independent of the choice of P(m c , M\m c , M ,T) 
and P(f,M\f,M,l). 

No other observed binary pulsar is known to have a 
ns companion. Both / and M have also been deter- 
mined for PSRs 1855+09 and 1802-07 [10,2]; however, 
in these systems the companion is thought to be a white 
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dwarf. The constraints these observations place on the ns 
mass depend on the uncertain distribution of white dwarf 
masses in pulsar-white dwarf binaries; consequently, I do 
not consider these systems in making my estimates. 

In addition to ns-ns binary pulsars, where the masses 
are determined through observations of the pulsar 
Doppler velocity curve with precision sufficient to observe 
relativistic effects, ns masses have also been determined 
in several X-ray binary systems [11,12]. These masses 
are estimated from Keplerian-order observations of both 
the pulsar and the companion Doppler velocity curves, 
and an estimate of the orbital inclination angle based on 
the eclipse duration. The derived masses are more uncer- 
tain than for any of the ns-ns binaries considered here; 
additionally, X-ray binary neutron stars may constitute 
a sub-population with a different mass distribution than 
ns-ns binaries. In order to maintain a homogeneous sam- 
ple of neutron stars as well as to avoid the complications 
inherent in the modeling of the inclination angle I have 
excluded the X-ray binary observations from the sample 
considered here (however, see sec. d). 

c. Probability and statistics of neutron star masses 
Could we examine all neutron stars, we would know the 
exact distribution of their masses. Lacking complete 
knowledge, we can examine an incomplete sample and 
a family of hypothetical distributions and calculate the 
conditional probability that the sample is drawn from a 
member of the family. As the sample size decreases or the 
observations become less precise, we are less able to dis- 
criminate between distributions that differ only slightly. 
For small or imprecise samples, only the grossest features 
of the distribution can be characterized. From observa- 
tions of four ns-ns binaries all we can realistically hope to 
learn of the ns mass distribution is its mean and extent. 
Here I consider a simple family of distributions that cap- 
tures these gross features. As more observations become 
available, our understanding can be refined by refining 
the family of distributions. 

Thus, suppose that ns masses are uniformly distributed 
between upper and lower bounds m u and m; , so that 
P(m\mi,m u ,T) = 1/ (m u — mi) when m; < m < m u , 
and otherwise. In this Letter I find P(mi,m u \{g},2), 
where {g} represents the ns-ns binary observations from 
which the masses are determined. 

Using Bayes Law of conditional probabilities, 
P(mi,m u \{g},X) may be factorized: 

P { mi,m u \{g},2)= ^M^j^P(mi,m u \2)^ (2) 

The probability density P(mi, m u \2) represents our prior 
knowledge of the upper and lower bounds m u and m; . 
Impose here only the most conservative theoretical con- 
straints: causality and general relativity together pro- 
vide that neutron stars cannot be more massive than 
M u ~ 3 Mq [13,14], and our understanding of the sub- 
nuclear density equation of state provides that they are 



not less massive than M; ~ 0.1 Mq [15]. As a result, 
M u > m u > m; > Mi . Theory providing no further clear 
guidance, assume that within these constraints all pairs 
(m;, m u ) are equally likely; thus 

P{m u m u \X) = 2/{M u -M l f (3) 

for Mi < m; < m u < M u and otherwise. 

For a given set of observations {g}, P({g}\2) is a con- 
stant whose value must be such that 

pM u i>M u 

1= / dm, / dm u P(m,,m u \{g},T). (4) 

J Mi J mi 

Thus, once P({g}\mi, m u ,T) is known, P(mi, m u \{g},2) 
is determined through this normalization integral and 
there is no need to find P({g}\2) separately. 

To evaluate P({g}\mi,m u ,2), note that the observa- 
tions of each binary pulsar system are independent; con- 
sequently, 

P({g}\m h m u ,2) = Y\_P(g n \mi,m u ,2) (5) 

n 

where g n represents the observation of system n. The 
form of P(g n \mi,m u ,2) is determined by the charac- 
ter of the observation g n . For PSRs 1913+16 and 
1534+12, we need the joint probability distribution 
P(m c , M\mi,m u ,2) of the observed companion mass m c 
and total mass M for fixed m; and m u . Using Bayes 
law, this distribution can be expressed as an integral over 
P(m c ,M\m c ,M,2) (cf. eq. 1): 

P(m c , M\m h m u ,2) = 

// (Klt( ,MP (mt ,M|,7 lt ,M.I)P( ftt .M|,,,,,,,,,,l) 

(6) 

where P(m c , M\mi, m u , 2) = (m u — m{) 2 . 

For PSRs 2127+11C and 2303+46, we need the distri- 
bution P(f, M\mi, m u ,2) of the observed mass function 
/ and total mass M for fixed m; and m u . This distribu- 
tion can be expressed as an integral over P(f, M\f, M ,2) 
(cf. eq. 1): 

P(f,M\m h m u ,2) = 

J J dfdM P(f,M\f,M,2)P(f,M\m h m u ,2) (7) 

To evaluate P(f ', M\mi,m u ,2), write it as 

P(f,M\mi,m u ,2) = 

P(f\M,m,,m u ,2)P(M\m l ,m u ,2). (8) 

The two probability densities on the right-hand side can 
be calculated separately: 
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P(M\m,,m u ,T) = 

max , min ( M — mi, m u ) — max ( M — m u , mi 



(m u - miY 



P(f\M,mi,m u ,X) 



(sin 1 x\ — sin 1 xq) [M/f 



2/3 



1 

^ min [ M — mi, m u ) — max ( M — m u , mi 
where 



(9a) 



(9b) 



i n = max 



x-, = max 



0,1 



0,1 



2/3 



min | m u , max \ M — m u , mi 

\ 2/3 

max I mi , min ( M — mi, m u 



(9c) 



(9d) 



In the event M is unknown but we believe the system 
consists of two neutron stars, we can write 

P(f\mi,m u ,l) = J dM P(f,M\mi,m u ,l) (10) 

Either of expressions (8) or (10) may be useful in in- 
terpreting observations of binary pulsars systems in the 
context of a model for the pulsar mass distribution. 

d. Discussion Figure 1 shows contours of constant 
P(mi,m u \{g},X) containing 95% (solid line) and 68% 
(dotted line) of the total probability. Figures 2 and 3 
show the same contours if the survey is restricted to those 
systems where the individual masses are known (PSRs 
1913+16 and 1534+12) or to those where only the mass 
function and the total mass are known (PSRs 2127+11C 
and 2303+46). 

Summarizing figure 1, the 95% confidence intervals 
for the lower and upper bounds m; and m u when con- 
sidered separately, are 1.01 < m;/M < 1.34 and 
1.43 < m u /M.Q < 1.64, and the maximum likelihood 
values of m; and m u are 1.29 Mq and 1.45 Mq. The 95% 
confidence intervals for the mass range A = m u — mi 
and the total mass M , when considered separately, are 
0.10 < A/M < 0.51 and 2.46 < M/M Q < 2.88, and the 
maximum likelihood values of A and M are 0.16 M and 
2.73M . 

Quite independently of the observations, numerical 
calculations of gravitational collapse supernovae have 
led Woosley and Weaver to "predict" ns masses in the 
range 1.15-2.0 Mq, with a preponderance between 1.3 
and 1.6 Mq [1]. This prediction is certainly consistent 
with the range derived here. 

While I have excluded X-ray binaries from the anal- 
ysis presented here (cf. sec. b), it is worth noting that 
the inferred mass of the Vela X-l ns is greater than 
1.50 Mq with 95% probability [11]. An upper limit m u 
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FIG. 1. Assuming ns masses are uniformly distributed be- 
tween mi and m u , observations of PSRs 1534+12, 1913+16, 
2127+11C, and 2303+46 determine the joint probability dis- 
tribution for mi and m u . Shown here are contours enclosing 
regions of 68% (dotted) and 95% (solid) of this distribution. 
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FIG. 2. As in figure 1, except that the contours are based 
on the constraints provided by observations of PSRs 1534+12 
and 1913+16. 
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FIG. 3. As in figure 1, except that the contours are 
based on the constraints provided by observations of PSRs 
2127+11C and 2303+46. 

greater than this is consistent with our results even at 
the 68% confidence level; however, m u greater than the 
most likely value of the Vela X-l ns mass estimated by 
[11] (1.85 Mq) lies outside the 95% confidence interval 
estimated here. 

LIGO will be extremely sensitive to the "chirp mass" 
[M = (mim2) 3/,5 (mi + raj) -1 ' 5 ] of an inspiraling binary 
system, but less sensitive to the individual masses [4,5]. 
In the same manner that our prior knowledge that ns 
masses could in no case be greater than 3 M or less 
than 0.1 M played a role in this analysis (cf. eq. 2, 3 
and intervening text), so our understanding of the range 
and distribution of ns masses gained from this exercise 
will play a role in determining the confidence intervals 
for measurements of M, mi, and rri2 using LIGO obser- 
vations. An accurate assessment of our prior knowledge 
is especially important in determining when a signal is 
sufficiently strong that it refines our understanding as 
opposed to affirming our existing prejudices. 

The ns sample used in this survey is highly selected: 
only components of ns-ns binary pulsar systems are in- 
cluded. In addition to the exclusion of X-ray binaries and 
isolated neutron stars, this means that half the sample 
are pulsars. It is believed that isolated, non-millisecond 
pulsars are a fraction 10 -4 of neutron stars, and that ns- 
ns binary pulsars number approximately 1/10 of these 
[16]. Since the mass distribution of the ns sub-population 
considered here may not be representative of neutron 
stars generally, application of these results to isolated 
neutron stars must be made cautiously. Without good 
estimates of the how the mass distributions of these dif- 



ferent populations may differ it is not possible to estimate 
the effects of these selection biases. Nevertheless, it is 
clear that the homogeneous subset of neutron stars con- 
sidered here has, with high confidence, a range of masses 
restricted in a way that our understanding of ns and bi- 
nary system formation and evolution do not, but even- 
tually must, confront. 
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